Simulations of cubic-tetragonal ferroelastics 
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We study domain patterns in cubic-tetragonal ferroelastics by solving numerically equations of mo- 
tion derived from a Landau model of the phase transition, including dissipative stresses. Our system 
sizes, of up to 256 3 points, are large enough to reveal many structures observed experimentally. Most 
patterns found at late stages in the relaxation are multiply banded; all three tetragonal variants 
appear, but inequivalently. Two of the variants form broad primary bands; the third intrudes into 
the others to form narrow secondary bands with the hosts. On colliding with walls between the 
primary variants, the third either terminates or forms a chevron. The multipy banded patterns, 
with the two domain sizes, the chevrons and the terminations, are seen in the microscopy of zirconia 
and other cubic-tetragonal ferroelastics. We examine also transient structures obtained much earlier 
in the relaxation; these show the above features and others also observed in experiment. 



PACS numbers: 81.30.Kf, 62.20.Dc 
I. INTRODUCTION 

Ferroelastics 1-3 are crystalline solids that undergo a 
shape-changing, structural phase transition with decreas- 
ing temperature T. The high-T or parent phase distorts 
spontaneously at the transition temperature T c to form 
one or more products or variants with identical energies 
but different orientations. In cubic-tetragonal (C-T) fer- 
roelastics for example, only one of three four-fold axes is 
retained below T c , giving three variants. 

Due to constraints imposed by neighbouring material, 
a single variant is found only rarely, perhaps only in very 
small grains; the displacement resulting from a strain of 
10~ 4 in a grain of size 10 //m (both typical values) is 
an order of magnitude too large to be accommodated 
by atomic rearrangements at grain boundaries. Multi- 
ple variants can arise also from independent nucleation 
events. An external stress can move the walls that sepa- 
rate the variants, so converting one variant to the othcr(s) 
at little cost in energy; the process has been observed 
for example by neutron diffraction of zirconia 4 . Related 
phenomena are strongly hysteretic stress-strain relations, 
shape- memory effects, etc 3 . 

Pockets of the parent phase can persist below T c , until 
the gain in condensation energy overcomes the cost of in- 
troducing domain walls. Compositional inhomogeneities 
combined with a strong dependence of T c on composition 
can also smear the transition, in cases to over 100 K; it 
is common to speak instead of a transformation. 

Domain patterns in ferroelastics have nothing in com- 
mon with those in conventional order-parameter sys- 
tems. Two-dimensional (2d) patterns in tetragonal- 
orthorhombic (T-O) ferroelastics like YBa2Cu30y_5 re- 
semble not at all those in 2d Ising models, though both 
have two variants; neither can 2d patterns in hexagonal- 
orthorhombic (H-O) ferroelastics be understood from 2d 
three-state Potts models, both with three. The differ- 
ence arises because (a) in order to maintain a coherent 



interface (no dislocations or disclinations), ferroelastic 
domain walls rotate the variants as well as join them, 
and (b) order-parameter strains alone are insufficient to 
understand the patterns. Disclinations, which have no 
counterpart in Ising and Potts systems, are generated 
however by wall collisions; they dominate patterns espe- 
cially in H-O-like systems (e.g. trigonal-monoclinic lead 
orthovanadate) . 

Much of ferroelastic theory follows Barsch and 
Krumhansl 5 in expanding the free-energy density in pow- 
ers of the strains and their derivatives. Analytical results 
are possible in a few cases 5 ' 6 , but structures formed by 
colliding walls require numerical effort for their under- 
standing. Such studies 7 " 9 in 2d gave static and transient 
domain patterns that reproduced nearly all aspects of 
those observed in H-O-like and T-0 ferroelastics. Ref. 
10 obtained similar results. 

An alternative approach, phase-field theory 11 , pre- 
dates Ref. 5 and has also been used extensively to un- 
derstand domain patterns in H-O-like materials 12 , C-T 
materials 13 ' 14 , etc. Ref. 15 discusses differences between 
the two approaches. And other lines have been pursued. 

In C-T ferroelastics, domain walls joining two tetrag- 
onal variants lie optimally in the cubic 110 planes 16 and 
are then twin walls. In the solutions 5,17 ' 18 for the twin 
wall, the strains depart only locally from their bulk val- 
ues, decaying exponentially with distance from the wall. 
Of course domain walls need not lie in the cubic 110 
planes; they are then not twin walls, and the departures 
decay algebraically 19 . 

Optical and electron microscopy of zirconia 20-22 , 
leucite 23 ' 3 , barium titanate 2 and other 24 C-T ferroelas- 
tics reveals a variety of patterns. The structures of 
small and large grains are respectively lamellar (with 
only two variants) and banded (with all three) ; the latter 
seems unique to C-T materials. The lamellar <-» banded 
transition is understood 2 ' 25 , though the analysis is ap- 
parently limited to variants oriented at n/2. Other as- 
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pects of the domain patterns have also been explained 
analytically 20,26,2 ^ 25 . The character of the patterns de- 
pends also on the thickness of the sample and on whether 
the surface examined is part of clamped specimen, or a 
free surface, or representative of the bulk 2,25 . The highly 
sensitive technique of birefringence imaging 27 promises 
to reveal further details of these patterns. 

The following presents results from simulating the time 
evolution of C-T materials. We obtain the equations 
of motion by expanding the free-energy density and the 
Rayleigh dissipation density to lowest possible order in 
the strains and their derivatives. We solve numerically 
for the displacement, using periodic boundary conditions 
and omitting the inertial term. Our late-time structures 
reproduce most features of the banded patterns found in 
C-T materials 2,3,20 " 24 . These structures were not found 
in the smaller systems used in the previous C-T simula- 
tions of Refs. 28 and 29; some were found in a phase-field 
study 14 . Our transient structures show in addition other 
wall configurations observed experimentally. 



II. LANDAU THEORY AND EQUATIONS OF 
MOTION OF ELASTICS 

The displacement u(r) of a material point is defined 
relative to its position r in the parent phase. The sym- 
metric strain tensor in this Lagrangian description is 
r]ij = \{uij + Uj t i + Uk,iUk,j), where u itj = dujdxj. We 
neglect the nonlinear term in r\ because it has no known 
qualitative effect on the domain patterns. The strains 
(all of which vanish in the parent phase) are defined by 
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these definitions differ slightly from those in Ref. 18. The 
deviatoric strains e 2 and e 3 form the two-component or- 
der parameter of the transition. The other strains e\ (the 
dilatational strain in the small-strain limit) and the shear 
strains e±, e$ and e$ are identically zero in the uniform 
product phase; they are required however to understand 
domain patterns, even for a single twin wall 5,17,18 . 

The six strains are obtained from the three components 
of the displacement u and so are not independent when 
they vary spatially; the second derivatives of the strains 
are linked by compatibility relations, necessary and suf- 
ficient conditions that the strains be derivable from u. 
We satisfy these relations implicitly by working with the 
components m. Refs. 30, 29, 10, which work directly with 
the strains, satisfy them explicitly by imposing them to 
obtain nonlocal relations between the order-parameter 



strains ei and e 3 ; the anisotropic, oscillatory nature of 
the kernels, obtained also in Ref. 15, provides much in- 
sight into domain structures and their relaxation. Nonlo- 
cal relations were developed much earlier 11 , though there 
appear to be differences 15 . 

In the Landau expansion of the free-energy density in 
the strains and their derivatives, the cubic symmetry of 
the parent phase permits three invariants to second order 
in the strains, e 2 , e 2 +e§ and e 2 +e|+e|. The correspond- 
ing stiffness coefficients A\, A 2 and A4 are linear com- 
binations of the Voigt coefficients. The order-parameter 
stiffness A2 softens with decreasing T as A2 = a(T — To). 
To describe the phase transition, one adds a term cubic 
in e2 and e 3 (this term breaks the rotational symmetry 
in (e 2 ,e3) space), and a quartic term for stability. The 
minimal density, that contains only essential terms, is 
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the last term gives the domain-wall energy (which pre- 
vents the system from dividing into arbitrarily small 
domains). We omit all unnecessary terms, namely 
higher-order terms and also some of the same order as 
those kept (a second invariant in the order-parameter 
derivatives 18,29 and other derivative invariants). 

The coefficients A 2 , B 2 and C 2 determine the transi- 
tion temperature T c and the spontaneous strain e 30 in 
the product phase. When A 2 > B 2 /4C 2 , the free en- 
ergy has only the cubic minimum at e 2 = e 3 = 0. For 
A 2 < B 2 /AC 2 , it has also three degenerate tetragonal 
minima located at 
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The phase transition, which is first-order, occurs when 
A 2 = |£> 2 /C 2 . The cubic phase is unstable for A 2 < 0. 
The symmetric stress tensor <Ty is defined by 

aij = ST/Srjij = GkSek/Sruj (12) 

with Qu = Sf/Sek; explicitly, 

Q 2 = (A 2 - D 2 W 2 ) e 2 + 2B 2 e 2 e 3 + C 2 e 2 (e 2 + e 2 ) (13) 

G 3 = (A 2 - D 2 V 2 ) e 3 + B 2 (e 2 . - e 2 ) + C 2 e 3 (e 2 + el) 

(14) 
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Stresses arise also from dissipative mechanisms. The 
same symmetry considerations as used for the free energy 
give the Rayleigh dissipative density as 

* = f ej + ^(e| + e\) + ^(e| + ef. + e%) (21) 

to lowest order in the time derivatives ej of the strains. 
The dissipative stresses a'^ are found from 

4 - 5H/5fHj = S' k 6e k /5fHj (22) 

where ^ = <5ft/<5e fc : # = A[e u ■■■,G' 6 = A' 4 e 6 . 

Our interest is in static states and in states where walls 
move slowly (rather than in effects associated with mo- 
tion at or near the sound velocity) and so we assume 
isothermal conditions. The equations of motion follow 
from Newton's second law 

fi = pu t = a'^ + (Tijj . (23) 

In the overdamped limit, the inertial term p%n is dropped 
and Eq.(23) simplifies to a' i3 ■ • = In terms of the 

displacement, this is 

[A'df + B' (a| + flf)] ui + C (8^2 + d!d 3 u 3 ) 
= -(A- \D 2 V 2 ) df Ul - B (dl + flf ) ui 

- (C + \D 2 V 2 ) (ftftua + didsua) + R? L (24) 

for i = 1, with obvious forms for i = 2, 3. The coefficients 
are A' = A[ + \A' 2 , B' = \A' A and C = A\ - \A' 2 + \A' 4 ; 
the definitions for A, B and C are obtained by dropping 



the primes. The nonlinear terms 

R^ = -8,(1-0^ + ^0^) (25) 

R^ = -d 2 (-\G^ L + ^ L ) (26) 

< L - -ds (-^ L ) (27) 



on the right-hand sides involve the nonlinear parts of G2 
and Q 3 , namely the terms with coefficients B 2 and C 2 in 
Eqs.(13) and (14). 

The matrix on the left-hand side of Eq.(24) must be 
invertible. Special handling is required when the strains 
are constant (that is 9i = $2 = $3 = 0) since both sides 
are then zero. Examining the full equation of motion 



(23), we see that the proper way to account for the con- 
stant strains is to leave them constant at all times. The 
interpretation is physical: a piece of strained material 
does not move unless there is a differential strain. An- 
other case of interest is A[ = A' 4 = 0, which may be 
a reasonable choice given that it is the order parame- 
ter, and not the other strain components, which changes 
most quickly in time. Then the matrix is not invertible 
for di — or d 2 = or <9 3 = 0, but again Eq.(23) tells us 
how to handle this case. 

The equations of motion (24) can be solved under a 
variety of boundary conditions. Wishing to examine do- 
main patterns, we used periodic boundary conditions 
(which allow no length change in any direction) in or- 
der to force domain walls into the low-T phase; all three 
variants are required since two variants can satisfy the 
constraint in only two directions. A second important 
consideration is that these conditions allow use of the 
fast Fourier transform, which is much faster than real- 
space methods. We should however voice our concern 
that these conditions may lead to spurious correlations 
between relaxation events at large relative distances. Of 
other choices, clamped conditions (u = on and out- 
side the boundaries, as in Ref. 8) would also force do- 
main walls into the low-T phase, but are less attractive 
because they usually give complex structures near the 
edges. Open conditions are not useful for our purpose, 
for the system would go to a single variant for almost 
any initial state. Yet another possibility corresponds to 
applied stresses at the boundaries. 

In Fourier space, Eqs. (24) are identical to the 
equations of time-dependent Ginzburg-Landau (TDGL) 
theory 30 ' 29 ' 10 . In real space, Eqs. (24) contain extra 
space derivatives on both sides relative to TDGL theory 
and so appear more general (they can be applied for ex- 
ample to systems with imposed strains). The neglect of 
the inertial terms both here and in Refs. 30, 29, 10 is 
problematic for the relaxation. 

III. DOMAIN PATTERNS 

We solved Eqs. (24) numerically, as described in the 
Appendix, starting usually from random displacements. 
We present fully converged 31 results for 3d grids of 
128 3 points and quasi-static 31 results for 256 3 points. 
We present also transient structures in systems of 256 3 
points. Systems of these sizes reveal features not seen in 
previous studies, which used at most 64 3 points. 

The nature of the patterns depends in part on the stiff- 
ness coefficients A\ and A4. The dilatational and shear 
energies arc minimized when the walls lie in cubic 110 
planes, and so the parameters A\ and A4 control the 
energy cost incurred when the walls depart from their 
optimal orientations. In stiff systems, with large A\ and 
At, the domain walls must lie close to the 110 planes, 
whereas in soft systems they can depart from these opti- 
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mal orientations at small cost in energy. We have how- 
ever no quantitative way to distinguish stiff from soft 
systems; comparing A\ and A\ with the order-parameter 
stiffnesses 7 found from the curvatures of the free energy 
about the tetragonal minima is not an effective means. 

A. Late-time structures 

We made extensive studies of late-time domain struc- 
tures at A 2 — — 2 and A 2 = —20; the spontaneous strain 
e3o at these temperatures is respectively twice and four 
times the value at T c . We present results at only the 
lower (latter) value, where the order parameters inside 
the domains are more well developed. 

Most of w 200 simulations gave multiply banded or 
herringbone structures like those shown in Figure l 32 . 
These patterns, which seem at first glance to reflect more 
the periodic boundary conditions than any physics, are 
in fact found in the microstructure of polydomain zirco- 
nia (examples are Figures 4(a), 5(a), 6 and 7 of Ref. 20, 
Figure 3 of Ref. 21, and Figures 1 and 2 of Ref. 22) and 
other C-T materials 24 . Similar patterns appear also in 
the well known ferroelectric BaTiC-3, specifically Figures 
2(a) and 8(b) of Ref. 2; they should appear in other elas- 
tic/electric and elastic/magnetic ferroics provided that 
the elastic energy dominates the electric and magnetic 
energies, as it does 2 in BaTiC^. Banded structures were 
found also in a phase-field study 14 . 

All three variants appear in Figure 1, but not equiva- 
lently. The structures consist of two primary bands, here 
red and green; the width of the primaries is determined 
by the system size in Figure 1 and, one assumes, by the 
grain size in experiment. Each primary is penetrated by 
the third variant, here blue; neither primary contains do- 
mains of the other. Within each primary, the host and 
the third variant form secondary bands; the ratio of the 
width of the host variant to the width of the third variant 
in the secondary bands is ideally 2:1 so that the three 
variants appear with equal volume fractions 20,26 . The 
same ratio was found in Refs. 2, 25, which found also the 
optimal value of the period of the primary bands to that 
of the secondary bands. 

Figure 1 (a) shows a fully converged 31 structure; it has 
the lowest energy of seven states found in 24 quenches 
with the same parameters. The front face shows chevron 
(or herringbone) structures; the blue domains are contin- 
uous across the red-green boundaries, where they bend 
through 90° and are slightly distorted as well. 

Figure 1(b) shows another fully converged structure 
obtained with the same parameters as part (a); it has 
a higher energy (the third lowest of the seven states) 
and so is metastable. Some of the blue variants form 
chevrons, as in part (a), but some terminate at the red- 
green boundaries; the walls occasionally deviate from the 
110 planes. Presumably the higher energy relative to 
part (a) results in part from the terminations and the 



deviations; both are seen experimentally, for example in 
Figures 4(a), 5(a) and 6 of Ref. 20. 

Figure 1(c), for a 256 3 system with stiffer parameters, 
shows a quasi-static 31 configuration with a mixture of 
continuing and terminating variants. The system is large 
enough to show the secondary banding clearly, but it is 
too small and has too many imperfections to display well 
the 2:1 ratio discussed above. Of the two other simula- 
tions performed with the same parameters, one gave no 
terminations and not surprisingly a lower energy, and the 
third gave more terminations and a larger energy. 

In addition to these banded structures, we found lamel- 
lar structures in smaller systems (64 3 ), particularly for 
stiffer parameters; the agreement with the lamellar <-> 
banded transition analysed in Ref. 2, 25 is however 
largely illusory, for the order parameters cannot approach 
their optimal values due to the length constraint in the 
third direction. We found also intermediate structures 
in which the narrow variant appears in only one of the 
two primary bands. Finally, a few systems gave very dif- 
ferent tweed-like or basket-weave structures, of all three 
variants, that seem not to be observed in experiment. 

B. Transient structures 

In all our late-time banded structures (not just those 
of Figure 1), walls collide only at boundaries between the 
primary bands. Within each primary red (green) band, 

(a) the red-blue (green-blue) walls adopt only one of the 
two possible orthogonal orientations 16 , and 

(b) the other primary, the green (red) variant, is absent. 
Walls colliding within the primary bands are however 

observed; examples are the A3 band in Figure 6 of Ref. 
20, Figure 7(b) of Ref. 20, Figure 3 of Ref. 21, and Figure 
11 of Ref. 22. Because we obtain only two primary bands, 
neither did we observe the A2 and A4 bands in Figure 6 
of Ref. 20; these contain the same two variants as the Ai 
band but with the orthogonal wall orientation. 

Perhaps our systems are too small to show these ef- 
fects, perhaps the experimental systems are incompletely 
relaxed; that different experimental conditions can give 
different patterns 2 may be relevant here. We have how- 
ever found some of these features in transient structures 
of a 256 3 system; 128 3 systems are too small to show 
interesting features clearly. 

The important aspects of Figures 2(a) and (b) are the 
following. 

1. Needle twins: The top face of part (a) shows a band 
of green needle twins in the red primary band, and also 
a band of blue needle twins in the green primary band; 
needles appear also elsewhere. In part (b), some needles 
have advanced and some have retracted. Needle twins 
are found for example in leucite (Figure 3 of Ref. 23). 

2. Collisions of identical variants: In the green pri- 
mary band on the top face, green/blue walls collide with 
green/blue walls of the other orientation, forming modu- 
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lated structures. Figure 11 of Ref. 22 shows a similar pat- 
tern in zirconia, "a rather exceptional case" , also formed 
by orthogonal colliding walls. As in Refs. 8, 9, we ascribe 
these modulations, and also the structures in Figure 3(b) 
of Ref. 23, to formation of wedge disclinations between 
two identical but differently rotated variants; the same 
explanation applies to tip splitting, in some cases. 

3. A split tip: Tip splitting seems to occur only rarely 
in C-T materials (relative to T-0 materials), presumably 
because of the extra freedom afforded by three variants; 
an example is Figure 3 of Ref. 23. We found only one split 
tip, an indistinct one at that, at the right side of the top 
face in part (a); this is of course a transient configuration. 
Split tips appear in the statics of some simulations 8 ' 14 , 
but only at the interface with parent material; they are 
more frequent in transient structures 9 . 

4. Collisions of different variants: In the top faces of both 
parts (a) and (b), green needle twins collide orthogonally 
with, or come near, blue variants; collisions occur also 
in the lower front face of part (a). No special features 
result from the collisions, which have not been noted in 
any experiment known to us. 

Figure 2(c) shows the same system at the later time 
t = 30. Many of the defects in part (b) have disap- 
peared in this quasi-static wall configuration; the walls 
are straighter, but in this relatively soft system they still 
bend where the third variant (here blue) terminates. 

Inspection of the dilatational and shear strains of late- 
time structures shows, not surprisingly, that their magni- 
tude is maximum in the wall-collision regions. We inves- 
tigated also the early stages of growth initiated by locally 
perturbing the parent phase at a temperature well below 
T c ; the growth occurs predominantly along spikes in the 
111 directions and planes in 110 directions, producing a 
noncompact object. 
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APPENDIX A 

We solved Eqs.(24) using periodic boundary conditions 
on the displacement u. Since only qualitative compari- 
son with experiment seems possible at present, and be- 
cause we wished to obtain results qualitatively applica- 
ble to many C-T materials, we scaled the energy, the 



strains and the length; this scaling requires neglect of 
the nonlinear term in the strain tensor r\. We chose the 
values B 2 = 3 x 10 3 and C 2 = 2 x 10 6 so that the tran- 
sition occurs at A 2 = 1, and the scaled strain at T c is 
e3o(T c ) = 10~ 3 , an arbitrary value. We chose D 2 = 1 to 
set the scale for the domain- wall width 18 . The number 
of parameters in the free-energy density is then reduced 
from six to three, the scaled stiffnesses A\ and A 4 and 
the scaled T-like variable A 2 . Results for other parameter 
values are easily obtained by scaling back to the original 
variables. 

The viscosity coefficients A'^ appearing in Eq. (24) are 
not known from microscopic theory. Neither we expect 
can they be determined from experiments such as ultra- 
sonic attenuation that operate on time scales very differ- 
ent from those governing domain- wall motion (as distinct 
from the "twin cry" in some materials). Lacking experi- 
ments that might determine the relevant coefficients, we 
chose the unit of time so that A' 2 = 1; lacking a reason 
to do otherwise, we chose A[ = A' 4 = 1 also. 

Each time step began with the Fourier coefficents 
Uj(k,t) at time t. The left-hand sides of Eqs. (24) are 
linear in the strains and so the space derivatives are ob- 
tained by multiplication in Fourier space; our approxima- 
tions for the derivatives are described below. The linear 
terms on the right-hand sides are found in the same way. 
To obtain the nonlinear terms Rf L in Fourier space, we 
formed the Fourier coefficients of the strains e 2 and e3, 
transformed them to find the strains in real space, found 
the nonlinear terms Q 2 L and Q^ L by multiplication (re- 
placing the strains point by point), transformed the two 
terms back to Fourier space, and multiplied to obtain the 
space derivatives in Fourier space. The solutions were ad- 
vanced in time by an Euler step (usually At = 4 x 10~ 3 ); 
solution of three linear algebraic equations then gives the 
three components Uj(k, t+At) in Fourier space. To mon- 
itor the convergence, we found the energy, the root-mean- 
square order-parameter strains and right-hand sides of 
Eq.(24), etc, every 10 or 20 time steps. 

The above computational scheme requires storage of 
five matrices, three for the Uj (k, t) and two for the strains 
e 2 and e$ (or Q 2 L and Q^ L ). The fast Fourier transforms 
were performed using the Numerical Recipes routine 
fourn 33 , which deals with complex matrices. The full 
executable file for a 256 3 system requires 1.35GB of stor- 
age. Savings of about two in storage and execution time 
would be obtained by use of routines for real matrices, at 
though some expense in coding and clarity. We have also 
obtained static structures by conjugate-gradient mini- 
mization of the energy 33 . The latter method is preferable 
in some respects to solving the equations of motion, and 
we have used it to verify the correctness of some of our 
64 3 and 128 3 results; it requires however roughly 7 times 
more storage, well in excess of that available to us for 
256 3 systems. 

The first and second derivatives were obtained from 
obvious generalizations of the Id finite-difference approx- 
imations 
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If the Fourier coefficients are denned by 
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with period L = Nh in each variable, then the first and 
second derivatives are approximated by 

|- - ^- sin(2vrift/L) [3 + 2 sin 2 (Trjft/L)] (A4) 



FIG. 1. Late-time domain patterns. Fhe faces of each dis- 
play cube are cubic 100 planes; the red, green and blue regions 
correspond to the three tetragonal variants; the domain walls 
(black) lie optimally in cubic 110 planes. Parts (a) and (b) 
show fully converged 31 patterns (both at t ~ 70) obtained in 
128 3 systems using different starting configurations but oth- 
erwise identical parameters (Ai = A4 — 100); pattern (b) is 
a metastable state with higher energy than pattern (a). Part 
(c) shows a quasi-static 31 configuration in a 256 3 system at 
t — 140 with stiffer parameters (Ai = A4 = 500). The tem- 
perature parameter is A2 = —20 for all three parts. 



FIG. 2. Snapshots of a 256 3 system as it relaxes from ran- 
dom initial displacements. The faces are cubic 100 planes. 
The temperature parameter is A2 = —20 and the stiffnesses 
are Ai = A4 = 100. Parts (a) and (b) are transient structures 
at times t = 2.2 and 3.6 respectively; part (c) is a quasi-static 
pattern at t = 30. 



dx 2 



3h 2 



sin 2 (njh/L) [3 + sin 2 (njh/L)] (A5) 



in Fourier space. In obtaining second-derivative terms 
like df in Eq.(24), it is important to use Eq.(A5) rather 
than Eq.(A4) twice, for the latter vanishes at j = N/2. 

The space step size ft must be chosen as a reasonable 
compromise between the conflicting demands of large 
physical size L = Nh on the one hand and accuracy on 
the other. Our values ft = 0.5 at A 2 = —2 and ft = 0.25 
at A2 = —20 were established as follows. 

We first performed 48 quenches at A2 = —20 with 
Ai = A4 = 100 on systems of identical linear size 
L = Nh; 24 of these quenches used (A, h) = (128, 0.125) 
and 24 used (N, h) = (64, 0.25). The larger step size gave 
4 states, each of which was clearly identifed with a state 
found for the smaller ft. The energies of the 4 states com- 
mon to the two sets of quenches agreed to better than 1 
part in 2000 (relative to the uniform product phase) and 
the relative frequencies of occurrence were comparable. 
The smaller step size gave however 3 additional states, 
each once. Since ft = 0.25 is satisfactory at A2 = —20 
and since the variational wall width 18 scales as l/e3o, 
one expects ft = 0.5 to be satisfactory at A 2 = —2 where 
is half the value at A2 = —20. Less extensive tests 
carried out at A2 = —2 with A\ = A4 = 100 gave compa- 
rable results for (A, ft) = (128,0.25) as against (64,0.5). 
Similar tests with A\ = A4 = 1000 at both values of A 2 
gave the same conclusions. 

In passing, we remark that the number of metastable 
states found in the quenches of the previous paragraph 
and in the quenches used for parts (a) and (b) of Figure 1 
is smaller than found in a T-0 study 8 of clamped systems 
of comparable linear size; the difference is due to the 
different boundary conditions. 
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